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Abstract 

Authors' note: We rephrased the abstract. 

We employ KAM theory to rigorously investigate quasiperiodic dynamics in 
cigar-shaped Bose-Einstein condensates (BEC) in periodic lattices and super- 
lattices. Toward this end, we apply a coherent structure ansatz to the Gross- 
Pitaevskii equation to obtain a parametrically forced Duffing equation describ- 
ing the spatial dynamics of the condensate. For shallow-well, intermediate- 
well, and deep-well potentials, we find KAM tori and Aubry-Mather sets to 
prove that one obtains mostly quasiperiodic dynamics for condensate wave 
functions of sufficiently large amplitude, where the minimal amplitude de- 
pends on the experimentally adjustable BEC parameters. We show that this 
threshold scales with the square root of the inverse of the two-body scattering 
length, whereas the rotation number of tori above this threshold is propor- 
tional to the amplitude. As a consequence, one obtains the same dynamical 
picture for lattices of all depths, as an increase in depth essentially only affects 
scaling in phase space. Our approach is applicable to periodic superlattices 
with an arbitrary number of rationally dependent wave numbers. 

MSC: 37J40, 70H99, 37N20 

PACS: 05.45.-a, 03.75.Lm, 05.30. Jp, 05.45.Ac, 03.75.Nt 

Keywords: Hamiltonian dynamics, Bose-Einstein conden- 
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1 Introduction 

Bose-Einstein condensates (BECs) have generated considerable excitement in the 
physics community both because their study allows one to explore new regimes of 
fundamental physics and because of their eventual engineering applications. They 
constitute a macroscopic quantum phenomenon, and their analysis has already lead 
to an increased understanding of phenomena such as superfluidity and supercon- 
ductivity. Of particular interest are BECs in optical lattices (periodic potentials), 
which have already been used to study Josephson effects 0, squeezed states [49J, 
Landau-Zener tunneling and Bloch oscillations |4*3*j . and the transition between su- 
perfluidity and Mott insulation [HU HS1 121] • With each lattice site occupied by one 
alkali atom in its ground state, BECs in periodic potentials also show promise as 
registers for quantum computers |HZl EH] • 
Authors' note: We rephrased the next paragraph. 

In the present paper, we generalize recent work on near-autonomous dynamics in 
BECs to study quasiperiodic behavior in BECs in periodic lattices, which can 

have shallow, intermediate, or deep wells. We present our methodology and results 
in section 11.11 In section we discuss the physics of BECs and use a coherent 
structure ansatz to derive a parametrically forced Duffing oscillator describing the 
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spatial dynamics of the condensate. The quasiperiodic dynamics of parametrically 
forced Duffing oscillators is rigorously investigated in sectional Sections 13 . II and I3~2l 
then describe applications to BECs in periodic lattices and periodic superlattices, 
respectively. The KAM theorem used in this analysis is proven in section 0] Finally, 
we summarize and discuss our results in section El 

1.1 Methodology and results 

Authors' note: We rephrased the paragraph below. 

The spatial dynamics of standing waves in BECs in periodic optical lattices can 
be described by a parametrically forced Duffing equation, where the periodic forcing 
is given by an external potential due to the lattice [TU1 E21 E3| • This gives a 1~ degree 
of freedom Hamiltonian system. We use KAM theory and Aubry-Mather theory to 
study its quasiperiodic dynamics. Previous KAM studies in BECs took a heuristic 
approach and considered only near-autonomous situations (SHI E2] ■ The approach of 
this paper, however, is especially versatile in that shallow, deep, and intermediate 
lattice wells can all be considered using the same mathematical framework. That is, 
we consider the near-autonomous and far-from-autonomous settings simultaneously. 

Theorem El proves that for any (analytic) external periodic potential, any neg- 
ative two-body scattering length, and any chemical potential, one obtains mostly 
2-quasiperiodic dynamics for condensate wave functions of sufficiently large ampli- 
tude, where the minimal amplitude depends on the experimentally adjustable BEC 
parameters. In particular, the threshold amplitude is proportional to the reciprocal 
of the square root of the scattering length. Any 2-quasiperiodic wave function above 
the threshold has one fixed frequency and one proportional to its amplitude. We 
also demonstrate numerically that one obtains the same dynamical picture for lat- 
tices of all depths, as an increase in lattice amplitude essentially only affects scaling 
in phase space. These numerical results support the theoretically predicted scaling 
of the threshold amplitude. Our theorem applies to periodic superlattices with an 
arbitrary number of rationally dependent wave numbers. 

The system we investigate is given by a Hamiltonian of the form H(R,S,^) = 
\S 2 +U £), where \R\ is the amplitude of the wave function, R and S are conjugate 
variables, £ G R/Z, and £' = 1. The function U is a polynomial in R and 1-periodic 
in £. We also consider its Poincare map, defined to be the return map on the section 

£ = o. 

One can show that such systems have invariant tori sufficiently far from the origin 
(R, S) = (0, 0) even when they are far from autonomous (that is, even when H is 
a large perturbation of an £-independent system). In the present case, this implies 
the existence of invariant tori for any optical lattice depth. The key condition is 
that U(R,£,)/R 2 — > +oo as R — > ±oo, as this guarantees that the set of frequencies 
corresponding to rotation around the origin is unbounded; indeed, the frequency 
goes to infinity with the distance to the origin. The first result of this type can be 
found in the equivalent context of adiabatic theory Another strand goes back 
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(a) (b) (c) 




Figure 1: Phase portraits of the Poincare map P for the example of section I3TT] with 
horizontal axis R and vertical axis S. In this case, the Hamiltonian is of the form 
H = \S 2 — |[l + Vi cos(^)]-R 2 + ^i? 4 . Observe the invariant curves for large R in these 
figures. Generically, the invariant manifolds of the central saddle point intersect 
transversally, creating a homoclinic tangle and thereby implying the existence of 
horseshoes of measure zero. By conjecture, the closure of these horseshoes is a 
"chaotic sea" of positive measure, corresponding to the dots in the figures. As the 
size of the perturbation V\ is increased to +00, the size of each of the remaining 
integrable islands vanishes, but their total measure remains 0(1). The values of V± 
are (a) V l = 0.1, (b) V 1 = 0.5, and (c) V l = 1. 

to the question of boundedness of solutions (which is implied by the existence of 
invariant tori in this low- dimensional situation) These qualitative results have 
also been extended to more general mathematical settings ffi H l4"5 ] 147) ] IHH ] Ifi?) ] \'M \ \'<V2\ . 
Authors' note: We added Levi91 to the above list of references. 
Authors' note: The paragraph below was rephrased. 

A typical system of this type exhibits a phase space divided into two clearly 
distinct regions. See figures E El which show numerical experiments in the present 
setting. One region, bounded away from the origin, consists largely of invariant tori, 
with small layers of chaotic dynamics between them. Indeed, one can show that the 
measure of these layers vanishes exponentially fast as the distance to the origin goes 
to infinity |37| 15]. The other region has mainly chaotic dynamics. In this paper, 
we employ a quantitative existence result ^| to obtain bounds on the location and 
frequencies of invariant tori. 

2 Physical Background 

At low temperatures, particles in a dilute boson gas can reside in the same quan- 
tum (ground) state, forming a Bose-Einstein condensate (BEC) [T7]. This was first 
observed experimentally in 1995 with vapors of rubidium and sodium QB1- 111 
these experiments, atoms were confined in magnetic traps, evaporatively cooled to 
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(d) 



(f) 



Figure 2: Continuation of figure H Phase portraits of the Poincare map P at (d) 
Vx = 10, (e) Vi = 25, and (f) V x = 100. 



(g) 



(h) 



(i) 



Figure 3: Continuation of figure 121 Phase portraits of the Poincare map P at 
Vy = 200, (h) Vi = 1000, and (i) V x = 10000. 
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temperatures of a few hundred nanokelvin, left to expand by switching off the con- 
fining trap, and subsequently imaged with optical methods [17 \. A sharp peak in 
the velocity distribution was observed below a critical temperature, indicating that 
condensation had occurred. BECs are inhomogeneous, allowing condensation to be 
observed in both momentum and position space. The number of condensed atoms 
N ranges from several thousand to tens of millions. 

A BEC has two characteristic length scales: the harmonic oscillator length 
dho = ^/h/ [fTiLUho] (which is about a few microns), where ujho = (^x^y^z) 1 ^ 3 
is the geometric mean of the trapping frequencies, and the mean healing length 
X = 1 1 v/87r|aj^ (which is also about a few microns) , where n is the mean density 
and a, the (two-body) s-wave scattering length, is determined by the atomic species 
of the condensate. Interactions between atoms are repulsive when a > and attrac- 
tive when a < 0. For a dilute ideal gas, a « 0. The length scales in BECs should 
be contrasted with those in systems like superfluid helium, in which the effects of 
inhomogeneity occur on a microscopic scale fixed by the interatomic distance |17j . 

When considering only two-body interactions, the BEC wave function ("order 
parameter") \I/(r, t) satisfies the Gross-Pitaevskii (GP) equation, 

ifMt = ( —2^- + 0ol* I 2 + V(f)J * , (1) 

where \I/ = ^(r, t) is the condensate wave function, V(f) is the external potential, 
and the effective interaction constant is go = [A7ch 2 a/m}[l + 0(( 2 )}, where ( = 
a/I^PIciI 3 is the dilute-gas parameter jTHEHI- A BEC is modeled in the quasi-one- 
dimensional (quasi-lD) regime when its transverse dimensions are on the order of 
its healing length and its longitudinal dimension is much larger than its transverse 
ones fOJEl- in the quasi-lD regime, one employs the ID limit of a 3D mean-field 
theory rather than a true ID mean-field theory, which would be appropriate were 
the transverse dimension on the order of the atomic interaction length or the atomic 
size. The resulting ID equation is [d ED] 

h 2 

itiipt = -tt-^ + + V(x)1> , (2) 

2m 

where ip, g, and V are, respectively, the rescaled ID wave function, interaction con- 
stant, and external trapping potential. The quantity \if)\ 2 gives the atomic number 
density. The self- interaction parameter g is tunable (even its sign), because the scat- 
tering length a can be adjusted using magnetic fields in the vicinity of a Feshbach 
resonance [2H EDI • 

Potentials V(x) of interest include harmonic traps, periodic lattices and super- 
lattices (i.e., optical lattices with two or more wave numbers), and periodically 
perturbed harmonic traps. The existence of quasi-lD cylindrical ("cigar-shaped") 
BECs motivates the study of periodic potentials without a confining trap along the 
dimension of the periodic lattice [SZj. Experimentalists use a weak harmonic trap 
on top of the periodic lattice or superlattice to prevent the particles from spilling 
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out. To achieve condensation, the lattice is typically turned on after the trap. If 
one wishes to include the trap in theoretical analyses, V(x) is modeled by 



where k% is the primary lattice wave number, k 2 > «i is the secondary lattice wave 
number, V\ and Vi are the associated lattice amplitudes, and Vh represents the 
magnitude of the harmonic trap. (Note that V\, V%, Vh, K\, and K2 can all be tuned 
experimentally.) When Vh <C Vi,V2, the potential is dominated by its periodic 
contributions for many periods. BECs in optical lattices with up to 200 wells have 
been created experimentally [oTIj . 

In this work, we let Vh — and focus on periodic lattices and superlattices. 
Spatially periodic potentials have been employed in numerous experimental studies 
of BECs (see, for example, Refs jTJ EE]) and have also been studied theoretically 
(see, for example, [TUl H2 EH ESI)- In recent experiments, BECs have also been 
loaded successfully into superlattices with k 2 = 3/«i [HI] • Additionally, over the past 
couple years, there has been an increasing number of theoretical studies on BECs 
in superlattices |22l ESI EH ESI IE] ■ 

Authors' note: The next paragraph has been rephrased. 

As mentioned in section UTTl and proven below, we obtain bounds on the location 
and frequencies of invariant tori. Our paper establishes quasiperiodic dynamics for 
sufficiently large amplitude. By a slight adaptation of the argument (see remark [TUl 
below), it can be shown that the quasiperiodic dynamics has large measure at large 
amplitude. It is generally conjectured that chaotic dynamics exist where invariant 
tori have been destroyed. In the present model, the presence of chaos (which was 
studied in BECs in, e.g., Refs. [12] HU E2 EH]) would reflect an irregular spatial 
profile R(x) (where |-R(x)| = \ip(x, t)\ is the amplitude of the BEC wave function), 
whereas invariant tori correspond to regular (i.e., quasiperiodic) spatial profiles. 

map: how is the phrasing above? 

Remark 1 : When the optical lattice has deep wells (large \V\\ or IV2I), one can 
also obtain an analytical description of BECs in terms of Wannier wave functions 
using the so-called "tight-binding approximation" 02]- I n this regime, the BEC 
dynamics is governed by a discrete nonlinear Schrodinger equation, which is derived 
by expanding the field operator in a Wannier basis of localized wave functions at 
each lattice site. 

Coherent structures solutions are described with the ansatz 



where R 6 R gives the amplitude dynamics of the wave function, 6 gives the phase 
dynamics, and the "chemical potential" /i, defined as the energy it takes to add one 
more particle to the system, is proportional to the number of atoms trapped in the 
condensate. When the (temporally periodic) coherent structure (Jl]) is also spatially 
periodic, it is called a modulated amplitude wave (MAW) |EH] . 




(3) 




(4) 
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Remark 2 : The present work is concerned with the spatial amplitude dynamics of 
solutions of the form (J3J). To ensure that such solutions are physically relevant, it is 
important to examine their stability with respect to the dynamics of the GP equation 
((2)). Bronski and coauthors jHllElHII were able to obtain rigorous stability results 
in some situations using elliptic-function potentials. (For repulsive condensates, for 
example, they used the ansatz (j3J) to construct elliptic-function solutions to (J2J), 
whose linear stability they proved for R(x) > by showing that they were ground 
states of the GP equation.) Their results can be applied to trigonometric potentials 
by taking the limit as the elliptic modulus approaches zero. More generally, one can 
address the stability of solutions of the form (J3J) through direct numerical simulations 
of (0) using such solutions as initial wave functions: -^(^(J) = _R(x). It was shown 
previously using both lattice and superlattice potentials V(x) that one can obtain 
numerically stable MAWs for @ with solutions of the form (j3J) with trivial phase 
9 = jSH Furthermore, the numerically stable "period-doubled" solutions 

(whose spatial periodicity is twice that of the optical lattice potential) constructed 
via subharmonic resonances using the ansatz (@J) with trivial phase (8 = 0) have 
very recently been observed experimentally 



Inserting (jlj) into the GP equation (j2J) and equating real and imaginary parts, 
one obtains 
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hnR(x) = R"(x) + 

2m 



?L[e\x)f + gR\x) + V(x) 



R(x) , (5) 



= ^- [26\x)R'(x) + e'\x)R(x)} , 
which gives the following nonlinear ordinary differential equation: 

The parameter c is defined via the relation 

which plays the role of conservation of "angular momentum," as discussed by Bronski 
and coauthors [TO] . Constant phase solutions (i.e., standing waves) constitute an 
important special case and satisfy c = 0. In the rest of the paper, we consider only 
standing waves, so that 

R" = ^- + -^-R? + - w V(x)R . (8) 

Remark 3 : When V(x) = 0, the dynamical system (jHJ) is the autonomous, 
integrable Duffing oscillator. Its qualitative dynamics in the physically relevant 
situation of bounded \R\ is illustrated in figure |U The methodology developed in 
the present paper can handle attractive BECs (g < 0) with either fi < or \i > 
but not repulsive BECs, as equation (jHJ) has unbounded solutions when g > 0. 
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Figure 4: Phase portraits of coherent structures in BECs with no external potential. 
The signs of [i and g determine the dynamics of (JHJ). (a) Repulsive BEC with 
/i > 0. Orbits inside the separatrix have bounded amplitude |-R(x)|. The period of 
such orbits increases as one approaches the separatrix. In this case, the dynamical 
system can be rescaled so that R" = —R + R 3 . (b) Attractive BEC with fi > 0. 
The dynamical system can be rescaled so that R" = —R — R 3 . (c) Attractive BEC 
with fi < 0. Here there are two separatrices, each of which encloses periodic orbits 
satisfying R ^ 0. The dynamical system can be rescaled so that R" = R — R 3 . 



3 Main Result 

This section states the main theorem of this paper, which concerns the existence 
of quasiperiodic dynamics in a class of systems of the form (JHJ), including the cases 
with V(x) given by periodic lattices and superlattices. Applications to these cases 
are subsequently given in two subsections. 

The equation of motion (jHJ) describes a l| degree of freedom Hamiltonian sys- 
tem whose forcing is periodic. We will allow arbitrary analytic periodic potential 
functions V. We first rescale the period of the forcing to 1. We introduce the phase 
variable £ = T~ x x (mod 1), where T > is the minimal period of V. Letting ' 
denote differentiation with respect to £, we define S = R', and 

Mt) = T 2 - %V(x)) , z 4 = -T 2 |§ , U(R,0 = z 2 (0R 2 + z,R\ 
This gives the suspended dynamical system 

R' = S, 

«' = 1, (9) 

with Hamiltonian 

H(R, S, = ^ 2 + U(R, . (10) 
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Authors' note: We changed the following remark. 

Remark 4 : We restrict to analytic systems for convenience only. It is not a neces- 
sary restriction, as our result is based on Herman's translated curve theorem [2*7| l28j. 
which requires only C 4 smoothness. Furthermore, finite-smoothness KAM results 
can be obtained from analytic ones by inverse approximation |58j . The restriction 
to analytic systems allows one to use Cauchy's integral theorem to replace C k norms 
by C° ones. 

Remark 5 : Because Z4 > 0, the function U is a well; that is, U — ► +00 as 
R — ► ±00 for all £. It is even in R, which has important ramifications for the sizes 
of the perturbations in our subsequent analysis. (See, for example, lemma Hoi below, 
where the leading term in the nonintegrable parts F\ and F2 would have been order 
one instead of going to zero as R — > +00 were it not for this symmetry.) 

We introduce action-angle coordinates as follows. (The details are in section 
H) Let H (R,S) = \S 2 + z 4 R 4 . For h > 0, define the action I = 1(h) to be 
the area in the (R, S^-plane enclosed by the curve H (R, S) = h. Let the angle 
= <p(R,S) G R/Z be such that the transformation (R, S) 1— > (0, /) is symplectic. 
This defines uniquely if we set 0(0, S) = for S > 0. 

In action-angle coordinates, the Hamiltonian takes the form K((j), I, £) = K (I) + 
Ki(<f>,I,£), where K (I) = H (R,S). We consider K as a perturbation of K . For 
any 1$ > 0, the unperturbed system Kq has an invariant torus I — Iq with frequency 
cj = Kq(Iq) in 0. We say that this frequency is of constant type with parameter 7 > 
if 

> 7 g~ 2 for all -GQ. (11) 
Q 

This is a special type of Diophantine condition [23 12H] ■ 

For a function / defined on a set V, we define ||/||x> = sup-p |/|. If / is vector- 
valued, then is the maximum of the norms of the components. For d > 0, let 
P(<0 = {£eC/Z:|Im(OI<2d}- 

With the additional notation 77 = I87, M = Z4I0, c = 126/25 = 5.04, and 
bi = J (l — u 4 ) 1 / 2 du = 0.874019..., the main result of this paper can now be 
stated as follows. 

Theorem 6 The Hamiltonian K has an invariant torus with frequency to = K' q (Iq) 
in if there exist v, j, d > 0, and 62 > with < v < | ■ 2 -7 / 3 and < 7 < || 



oj — 



P 
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such that 00 is of constant type with parameter 7 and the following conditions hold: 



A 



1 + £ 2 - 4/ V 3 >og(2), (12) 

1 < M, (13) 

47 

L < , (14) 

~ 200 ' v ; 

b 2 < 18[1 + 24&i(?7 + 2d)](l + L) 2/3 (l-L) 1/3 log(M), (15) 

2 < 5M 1/3 , (16) 

"* lllS » " (17) 



where 



log(M) ' 



" 1/3 -6- 4/3 (l - L)- 5 / 3 L 2 log(M),2M- 1 / 3 log(M) + — ^ 4/3 6 : 



B = 2- 13 / 3 3 3 7 8 c- 1 ^ 4/3 6 2 ^ 6 f - — tttt + ^^6- 4/3 6 2 ^^^) , 

1 2/ Vlog(M) 3 1 log (M)/ 

5 = -^^6 2 6r 5/3 (l + L)- 2 / 3 [l + 246 1 (7 7 + 2d)]- 3 . 

TTie torus Hes m the region given by \I — Z 1 < &2^o[log(z4/o)] _1 (i§^ + p)- 

Remark 7 : Conditions (|12|) - (|16|) are satisfied for M sufficiently large, while 
(fTTj) is a restriction on z 2 and hence on V. The theorem implies that H^Hixd) can be 
taken roughly proportional to z^ 3 I^ 3 . Equation (|2^j) below shows that z 4 1 ^I^ 3 is 
proportional to the maximal R coordinate -R max on the torus H Q (R, S) = K (I ), so 
1 1^2 1 |s(<2) is roughly proportional to z^R^^. This fits with the numerically computed 
phase portraits in figures d - El 

In terms of physical parameters, this implies that for an attractive BEC with 
given scattering length a < and chemical potential fi loaded into an arbitrary 
periodic lattice of amplitude ||V|| (with any number of wave numbers), the wave 
function's spatial component R(x) is quasiperiodic with two frequencies if its max- 
imum -R max is large enough, its frequencies satisfy the Diophantine condition pip, 
and the amplitude \ \V\ \ of the lattice potential is sufficiently small. The lower bound 
on -R max scales as k/ yja[, where k = 2ti/T is the lattice wave number (and we recall 
that a is the two-body scattering length), whereas the upper bound on ||V|| scales 
as ai? max . All frequency ratios that are algebraic numbers of index 2 satisfy the 
Diophantine condition (for some 7 > 0). 
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Remark 8 : The conditions of the theorem imply that M should be larger than 
roughly 10 6 . Indeed, from A > 8 • 2M~ 1 / 3 log(M) and BM 1 ' 3 > 2, it follows 
that A[l + 3 log(-B)/ log(M)] > 48M" 1/3 log(2). Hence condition JEJ) implies that 
M > 2 16 3 3 bf « 10 6 . 

As an illustration, we choose M and the parameters b 2 , 7, rf, v based on a numer- 
ical experiment where we evaluate the conditions of the theorem on a Cartesian grid 
of 360 x 180 x 120 x 120 x 20 points in the cube [10 6 , 10 18 ] x [lO^ 6 , 1] x [10~ 3 , 10] x 
[10- 3 , 10] x [10- 3 , § • 2~ 7 / 3 ] in (M, b 2 , 7, d, v) space, taking a logarithmic scale in the 
first four components and a linear scale in the last. For each M, we compute the 
largest value of the coefficient S over all grid points where all the conditions hold in 
order to obtain a good choice of parameter values. Figure El shows the largest S one 
can obtain for a range of M. 



logio ( S ) 




8 10 12 14 16 18 l g iQ ( M ) 

Figure 5: Log-log plot of M against the largest value of S found by a numerical 
computation. 



Corollary 9 If K has an invariant torus with frequency K' (Iq), then the return 
map of K on the surface of section £ = has an Aubry- Mather set with rotation 
number uo for any uo > K' (Iq). 

The corollary follows immediately from Aubry-Mather theory 01 0U1 IB] and the 
monotonicity of K' Q . Indeed, the Poincare map Pq corresponding to K is given by 

Po : + *£(/),/). 

From the monotonicity of K' Q , it follows that this is a twist map. The proof of the 
main result of Ref. shows that the Poincare map P corresponding to K also 
satisfies the twist condition. Alternatively, this also follows directly from the fact 
that the function U is superquadratic. 
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Authors' note: The paragraph that was here has been edited and changed 
into the following remark. 

Remark 10 : The translated curve theorem used (in (Ej) to prove theorem El shows 
persistence of invariant tori under relatively large perturbations (compared to other 
KAM or twist theorems), at the expense of excluding all but the "most quasiperi- 
odic" frequencies to. If the Diophantine condition is relaxed to include a nonzero 
measure of frequencies, then one can similarly show that the region outside the in- 
nermost invariant torus contains a set of invariant tori of positive measure, where 
the measure converges exponentially fast to full measure as / — > +00. The remain- 
ing dynamics in this region is presumably chaotic, but due to the low dimension of 
the system, there can not be any Arnol'd diffusion. Thus, it is at worst bounded 
chaos. Among many papers on this subject, see e.g. jHUHHllH] for further details. 

By contrast, the region on the inside seems to exhibit chaotic dynamics, as can 
be seen in figures \I\ El In the near-autonomous setting (i.e., for small-amplitude 
V), the figures show large domains of integrable dynamics in the interior region. A 
perturbation analysis for a similar system shows that, as the amplitude of V goes 
to infinity, the sizes of these islands vanish, but their number increases reciprocally, 
so that an integrable set of measure 0(1) remains |48j. 

In the next two sections, we consider potentials V(x) with one and two wave 
numbers. The former case describes BECs in periodic lattices jTUl H2 EE E2 155] . 
and the latter case, which is now experimentally accessible [5T], describes BECs in 
periodic superlattices [22l [HHl IMJ • 



3.1 Example 1: BECs in Periodic Lattices 

Optical lattice potentials are created experimentally as interference patterns of 
counter-propagating laser beams [T^j. In the periodic case, the external potential is 
typically taken to be sinusoidal, 

V(x) = -Vicos(Kx), (18) 

where k = 2ir/T is the lattice wave number. 
We write equation (jSJ) in the form 

R' = S, 

S' = —a.\R + a 3 R 3 + VlRcos(kx) , 

x' = l, (19) 

where a± oc u, oc g, and ' denotes differentiation with respect to x. The pa- 
rameters V\ , a± , CK3 G K and k > can all be adjusted experimentally. The phase 
variables are (R, S) G M 2 and x G M/(27tk _1 Z). The associated Poincare map P, 
which is defined to be the first return map on the section x — 0, corresponds to the 
flow over 2ttk~ 1 . 
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Equation (|19|) has two reversible symmetries: 

(R, S, x) I— > (-R, S, —x) , (P, S, x) I— > (P, — S, — x). 

The map P is thus reversible under reflection in both coordinate axes, 

R j oPoR j = p- 1 , j = 1,2, 

where Ri and R2 are the reflections with respect to the axes. 
Additionally, ([19)1 is invariant under two rescalings: 

(P, S, x; ai, a 3 , k, Vi) h- > (P, AS', A" 1 ^; A 2 «i, A 2 a 3 , Are, A 2 Vi) , (20) 
(P, S, x; oti, a 3 , k, Vi) 1— > (yuP, fiS, x; aii, fi~ 2 a3, k, Vi) . (21) 

The corresponding invariants for the Poincare map are obtained by dropping the x 
components. 

One can rescale ()19j) using these invariants to reduce to the cases where 
a\ = , ±1 and 0:3 = 0,±1. We look in detail at the case ai = —1, a 3 = — 1, 
corresponding to an attractive BEC with a negative chemical potential, where the 
underlying integrable system (with V\ = 0) is a "figure eight" (consisting of a central 
saddle point and two exterior centers, as shown in figure Eh). 

Because of the second rescaling (j21)l . the parameter a.3 simply measures the 
size of phase space. The first rescaling (|2T)|) shows that decreasing k increases the 
nonintegrable perturbation by a square law. Intuitively, a large lattice wave number 
k implies that the Poincare map corresponds to short- "time" integration in x and 
is thus "near" the vector field. With k rescaled to 1, as was done in stating and 
proving our main result, the perturbation V(x) in (fTH|) has unit period. 

Figures d - El show phase portraits of P at several amplitudes of the potential 
for (0:1,03, ft) = (—1,-1,1). One obtains qualitatively similar results for other 
values of the lattice depth V\ if the wave number k is rescaled, as indicated above. 
As remarked previously, the phase space in these phase portraits is divided into 
two clearly distinct regions: an outer one in which the dynamics consists in large 
measure of invariant circles and Cantor-like Aubry-Mather sets (that wind around 
the origin at large distance) and an inner one in which the dynamics is mostly 
chaotic. Our numerical simulations, which show the same scaling that our theoretical 
results indicate, suggest the presence of (parameter-dependent) integrable dynamics 
of positive but small measure inside the "chaotic sea." jlHj 

Remark 11 : A similar combination of islands of invariant tori within a chaotic 
sea occurs in the example of a parametrically forced planar pendulum The 
division of phase space into a mostly quasiperiodic and a mostly chaotic region is 
the typical behavior that one expects to observe in a large class of forced one dof 
Hamiltonian systems j!5j . 
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Applying theorem El for V(x) given by (jTSj) implies that the system ()19|) has an 
invariant torus with frequency vector (a;, 1) provided uo satisfies the conditions of 
theorem El The .R-amplitude is roughly equal to -R max , given by 

r) 'max — C)U 1 , j 

ny/-mg 

where g < for the present case of attractive BECs [see equation (|24|) and lemma 



3.2 Example 2: BECs in Periodic Superlattices 

Optical superlattices consist of small-scale lattices subjected to a long-scale periodic 
modulation. In recent experiments, BECs were created in superlattice potentials 
with a length scale (wave number) ratio of 1:3 |51j . However, theoretical research 
concerning BECs in superlattices has only begun to gain prevalence |2*2*1 I3b[ lo^j. 

To consider the case of (symmetric) periodic superlattices, we examine the po- 
tential 

V(x) = - [Vi cos(kxx) + V 2 cos(k 2 x)] , (22) 

where k 2 > K\ without loss of generality and k 2 /k\ G Q. The minimal period is 
T = 2tt/k, where k := gcd(fti ,K 2 ). 
Equation (jH} is then written 

— = S 

dx 

dS 

— — = —ctiR + a 3 R 3 + ViRcos(kix) + V 2 Rcos(k, 2 x) , 
dx 

where all the parameters are again experimentally adjustable. 

Applying theorem El for V(x) given by (|22|) with k 2 / ' K\ G Q (i.e., for periodic 
superlattices) implies that (|2^|) has an invariant torus with frequency vector (ui, 1) 
provided u satisfies the conditions of theorem El As in the regular lattice case, the 
i?-amplitude is roughly equal to -R max , where 

■KyJ-mg 

and we recall that 27t/k is the period of V(x) and g < for attractive BECs. 

Remark 12 : If k 2 /ki is not rational, then the potential V is not periodic. If 
Ki and k 2 satisfy a Diophantine condition, then one can prove the existence of 
invariant tori at large distance from the origin j22]- We conjecture that it is possible 
to quantify this existence result analogous to the periodic case jTHj. If /«i and k 2 
are not Diophantine — for example, if k 2 /k\ is a Liouville number — then one expects 
unbounded solutions and no invariant tori 
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4 Proof of the main result 



To prove Theorem El we first construct action- angle coordinates explicitly. We then 
employ the KAM theorem of Chow, et al. |T5J to complete the proof after a suitable 
transformation of the action variable. 

4.1 The action and Kq 

Define the action by 

1(h) = SdR = A ^2h - 2z 4 R 4 dR , 

JH =h JO 

where H (R,S) = ^S 2 + z 4 R 4 and -R max = Rmax(h) > is the solution of 
Ho(Rmax, 0) = h; that is, -R max = z^^h 1 / 4 . With the substitution u = R/R max 
and the relation h = -24-R max , the above integral reduces to 

1(h) = 4v / 2&izl /2 i?Lx = 4v / 26i^ 4 " 1/4 /i 3/4 . (24) 

Lemma 13 The unperturbed Hamiltonian in action-angle coordinates is given by 

K (I) = 2- 10 / 3 &- 4/3 zl /3 / 4 / 3 . 

The proof follows from the calculation above, noting that the unperturbed 
Hamiltonian K is the inverse of the function / [because H (R, S) = K Q (I) for 
I = I(H Q (R, S))]. The function I is invertible because 

g = Sv^^V 1 / 4 > 0. 
oh 

Remark 14 : One could also define the action to be the area enclosed by the 
curve Tj^ 2 + Z2R 2 + ZiR 4 = h, where z~2 is the average of z%, or even the area enclosed 
by H = h, where £ is considered as a parameter (that is, £' = in the unperturbed 
system). In the latter case, the action and angle will depend on £. Although 
these two approaches each leave a smaller term in the perturbation than our choice, 
and are therefore theoretically more pleasing, they lead to technical difficulties in 
the estimates we need to perform, as the expressions for the action and angle will 
involve elliptic integrals that depend on the phase variables. 

4.2 The angle 

In the upper half plane (S > 0), we define the angle by 

^^(lyTsoU*" (modl) ' 
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where S(h, R) = \/2h — 2z 4 i? 4 is the positive solution of H (R, S(h, R)) = h. Note 
that |^ = § H h S^ 1 dR. A similar definition holds in the lower half plane, where 

= \ — (ff ) _1 j$ S~ l dw (mod 1). Henceforth, we restrict to the upper half plane 
without loss of generality. 

Lemma 15 The following formulas hold: 

I rR/Ru 



cf>(h,R) = W J (1 - u*)- 1 ' 2 du , 



S(h,R) d 4(h,R) = 



y 246x i? max 

Defining r = r(/i, R) = R/ R mSuX (h), it follows that depends only on r: 

0(/i, J R) = $(r(h,R)), where 0(r) = / (1 - m 4 )- 1/2 dw . 

6oi Jo 

In particular, 0(1) = 1/4. 

Proof: Using the definition of and the substitution w = i? max M, we obtain 

which proves the first formula. Furthermore, 

S(h,R)^ = -(2h-2z 4 Rr /2 ^(l-rr i/2 RR^^ 



24&i 



24&i 



□ 



4.3 Localization and rescaling 

The system corresponding to K can be computed directly and is given by 

0' = ^(7) + / 1 (0,7,O, 
/' = / 2 (0,/,O, 
= 1, 

where h(cf>,I,0 = -Sg^ and ^(0,1,0 = 

17 



For a fixed Iq, we define a localization transformation (0,7, £) i— > (0, J, £) by 
I = I + ft(Io)J, with /3(/o) — ^-^oflogC^-^o)] -1 ) where 5 2 G K is positive and will be 
determined later. This takes the unperturbed torus / = Iq to J — and rescales the 
action variable, so that the components fi, / 2 of the perturbation are roughly the 
same size after rescaling and some other conditions are met (see remark l2~Tl below) . 
Although it is not a symplectic transformation, it nonetheless maps Hamiltonian 
systems to Hamiltonian systems. 

Define u = Kq(I ) and m = (3(Iq)Kq(I ), so that 

u = i. 2~« , (25) 

O 

m = \-^X' ! X^-. (26) 
9 log Wo) 

Lemma 16 In ((f), J, £) coordinates, the system is written 

fi^uj + mJ + gW + F^J,® , J' = F 2 (0,J,£) , f = 1, 
where, for some J*, J* G (0, J), 

= ^(/o) 2 ^"(/o + /9(/o)^)J 2 



J_ . 2 -4/3 fe -4/3 1/3.1/3 f x \ / 

27 61 " 4 /o I 1 + log Mo) ; I 



6 2 J 



27 V log (24/0)/ \logWo 

||(j) = p(i yKZ'(i +{3(i )j*)j 

1 , 1/3 7-1/3 / 7 T* \ -5,'" 

= L 2 -l/3 & -4/3 & ^4 j / x , °2 J \ 

27 1 2 log Wo) V log Wo)/ 

= |. 2 ^^ W v.( 1 + _W_) 

F 2 (0,J,O = /3(/ o )-V2(0,/o + /9(/o)J,O 

3 6, , 2 (e)r(l r) 6 2 + log Wo) J 



log (24/0 
1/3 



2/3 



4 J o 

— 1/3 

The lemma shows that the nonintegrable parts F\ and _F 2 are order I in 
leading term; this would have been order 1 had the original function U included a 
cubic term. 

Proof: Write Vt(J) = K' (I + (3(I )J), so that tu = £1(0) and m = Q'(0). By 
Taylor's theorem, 

g(J) = Q(J) -u-mJ = hi" (J*) J 2 , g'(J) = Of (J) - m = Q"(J*)J. 

Furthermore, (3(Iq)J' = I'. The expressions for g, g', Fi, and F 2 follow by direct 
calculation from (J2%|) and lemma ITol 
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□ 



4.4 Proof of Theorem [f3 



The rescaled system is a perturbation of 0' = u + mJ + g(J), J' = 0. The unper- 
turbed system has an invariant torus J = with ^-frequency u. Assume that u is 
of constant type with parameter 7 > 0. We now study the persistence of this torus 
under the perturbation F = (F\, F 2 ). 

For d > 0, 77 = I87, and p = (3m)~ 1 r], we define a (complex) neighborhood Vq 
of the unperturbed torus by 

V = V (r],p,d) = {(</>, J, £) E C/Z x C x R/Z : |lm(0)| < rj, \J\ < p} + 2d. 

The KAM theorem of Chow, et al. JH] now shows that the perturbed system has 
an invariant torus with frequency u satisfying \J\ < -^d + p provided 

dg 



d.J 



< m/4 . 

■Do 



(27) 
(28) 



The first of these conditions is a twist condition, and the second states that the 
perturbation is small enough, where c^o is defined as follows. Let c = 5.04, and 
||F||x> = max{||Fi||x) , H-F^Hdo}- Let Ly/ denote the Lambert W function (i.e., the 
inverse oiW ^ We w ). Then, 



^0 
b 

Q'O 

c 



a max 
2 + 3m 



log 2 
max < 12, 
3 + 12(77 + 2d) 



L w (b log 2) 

? 8 



108(7/ - 67) 4 
F||d + 2Cmax <j 1 



maxjmci, 2||F| \d } 
2\\F\ 



(29) 



Wo 



2d V ° I ' md 



■ max{m(2rf + p) + \ \g\\v + | \F\ \t> , 1} • 



Lemmas El ~~ 1201 show that the twist and smallness conditions follow from con- 
ditions (PU) - Cd of Theorem El 

Lemma 17 Condition \21\j follows from \14\j - 
Proof: For \ J\ < p + 2c?, it follows that 



b 2 J 
log(M) 



p + 2d 
log(M) 



L < 



47 
200 
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4 
2 



Figure 6: The set £, with labels for the four parts of its boundary (two semicircles 
and two line segments). 



Hence, for J* as in lemma ITol 

boJ* V 5/3 boJ 



1 + 



log(M) 



log(M) 



< 1 



47 
200 



-5/3 



47 3 
< - . 



200 8 



The desired result now follows from lemma fTol 



□ 



Lemma 18 Let r = R/R max , as before. Then, 

\\r\\ Vo < 1 + 24&i(t7 + 2d), 
||r(l-r 4 ) 1/2 |k < [1 + 2461(77 + 2d)] 3 . 

Proof: Without loss of generality, we restrict to (p satisfying Re(0) G [0, ■?] (mod 
1). Observe that |Im0| < rj + 2d in P - Thus, we can consider <p G T> = {<p G C/Z : 
Re(0) G [0, f](mod 1), \lmj>\ < rj + 2d}. 

Let £(e) = [0, 1] + e C C. It is obvious that r will be in £ for large enough 
e. Let e be the smallest number such that r = r(0) a £ = £(s) for all G T>. 
Hence, there exists 0o G P such that r(0o) is on the boundary of £. Furthermore, 
r(Re(0o)) € [0, !]■ Therefore, by the mean value theorem, 



e< |r(0o)-r(Re(0 o ))| < 



<9r 



(r? + 2d) . 



v 



By definition, §£(r) = (6bi)-\l - r 4 )" 1 / 2 . Thus, 

|^(0) = 66 1 (l-r 4 )V2. 

We will now show that ||1 — r 4 ||^ < (1 + e) A . By the maximum modulus theorem, 
||1 ~~ r4 |U * s attained on the boundary of £ , which looks like a stadium and consists 
of four parts (see figure |HJ). Part 1 is parametrized by r = eexp(i9), where 9 G 
[7r/2, 37r/2]. Hence, |1 — r 4 | < l+e 4 < (1 + e) 4 . Part 2 is parametrized by r = u — ei, 
with u G [0, 1]. Consequently, |1 — r 4 | < 1 — u A + 4w 3 e + 6u 2 e 2 + 4we 3 + e 4 < 
1 + Ae + 6e 2 + 4e 3 + e 4 < (1 + e) 4 . Part 3 is parametrized by r = 1 + e exp(i9), where 
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8 E [-tt/2, vr/2]. Hence, |1 - r 4 | < 4e + 6e 2 + 4e 3 + e A < (1 + e) 4 . Part 4 is similar 
to part 2. 

Thus, we have found that 

e < 6&i | |l - r 4 | |* /2 (77 + 2d) < 6^(77 + 2d)(l + ef . 

A straightforward calculation with this result then shows that e < 246j (77 + 2c?) , 
which implies the first statement of the lemma. To prove the second statement, we 
use the fact that ||1 — r 4 1 1 ^ < (1 + e) 2 < [1 + 246i (77 + 2d)} 2 , as already shown. 

□ 

Lemma 19 Assume that < TT3j) . (fT^l ), / TT3J/ 1 . and |77[) no/d. Tnen, 

9-4/3 

IMk < '—b-^il-LY^LHpl*, (30) 



Proof: The assumption on M implies that logM > 0. As in the proof of lemma 
IT71 ; b2 , J . r s < L < 1, and the result then follows from lemma ITol and lemma [PHI 

1 1 log(M) Vq — ' 1 1 1 1 

Condition (fTB]) implies that ||Fi||i5 < ||-p2||x> > an d (fT7|) bounds 1 1 ^2 1 1 x> - 

□ 



Lemma 20 Assume that v < | • 2 7 / 3 and 7 = n/18 < =| ; as m Theorem^ If 
conditions (Q3J) - fli7| ) ao/d ; taen 

cj n < - H ^— ^ M 1/3 . 

°" log(2) i3 + log(M)j 

Theorem |U] now follows immediately from this lemma and the expression for uo = 
K'(Io). 

Proof: By the previous lemma, 

C<-w a xfe+(— + V ) d6r 4/3 6 27 ^r + — &r 4/3 (l - LY^L'M^A . 
-3 \3 V9 y 1 log(M) 27 1 V ; ' J 

A straightforward calculation shows that v < |-2~ 7//3 implies 2||F||x> < md. Hence, 

3 + 12(n + 2d) n|F|| 

«o = ~ d \\W\W + 2C ) 



< 



{f M-V3 log( M) + + f ) d6 -/3 62+ 

- -6- 4/3 (l-L)- 5 / 3 L 2 log(M),2M- 1 / 3 /oa(M) + — cf6~ 4/3 6 : 1 



.4 



log(M) ' 
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From 7 = 77/I8 < 49/72, it follows that 12 < ^ • 7 8 (r] - 67)" 4 . We observe that 7 
does not appear in any other quantity besides b. The choice 7 = 77/18 minimizes b. 
We obtain 

7 8 • 18 2 ■ 3 4 (2 + 3m) , 

b = ; — - c rnd 

108 • 2 4 cr/ 6 

= 2~ 4 • 3 5 • 7 8 c- 1 7]-%2 + 3m)md 

= 2- 13 / 3 • 3 3 • 7 8 c~ l db- A/ % V ' e ( + —b~ 4 %J^-) M 1 / 3 

1 21 \\og(M) 3 1 2 \og\M)) 

= BM 1/3 . 

From b > 2, we see that Lw(b\og(2)) < log (6). The result follows. 

□ 

Remark 21 : The transformation 1 1— > J is chosen so that u and a>o have the same 
growth rate in I and 24 (i.e., in the leading term). Independent of the rescaling, 
io grows as I^ 3 . On the other hand, uj < a log (b)/ log (2) (if b > 2). Because b 
grows as some power of I , it follows that a must grow as I^ 3 / log Analyzing 
the dependence on z 4 then leads to the chosen rescaling. 



5 Summary and Conclusions 

In this paper, we discussed quasiperiodic dynamics in Bose-Einstein condensates 
(BECs) in periodic lattices and superlattices. The mean-field dynamics of a BEC 
is governed by the Gross-Pitaevskii (GP) equation which consists of a cubic 
nonlinear Schrodinger equation plus an external potential that takes into account 
the "trap" where the condensate resides. In this mean-field description, a given 
particle in the BEC is affected by the other particles only through average effects, 
and the non-local term in the original many-body Hamiltonian leads to the nonlinear 
term in the GP equation. One obtains a cubic nonlinearity if considering only 
two-body interactions. The GP equation, which is derived as a zero-temperature 
theory, provides a good description for BEC dynamics below the critical transition 
temperature at which the condensate forms [TT1 13*3*] . In "cigar-shaped" BECs, two 
dimensions are tightly confined, so one may further reduce ([*[| to ((2J), which has 
one spatial dimension [T7j. The amplitude dynamics of coherent structures of © 
with trivial phase (which describe standing waves) are governed by a forced Duffing 
equation given by (jSJ), or equivalently by (jHJ). The first equation arises directly from 
the physical setting, whereas ([""J) is the more convenient decription for mathematical 
analysis. We briefly describe the dynamics of this system. 

In the absence of the forcing (i.e., when the external potential V = 0), the 
system reduces to the autonomous Duffing oscillator and is integrable. Its dynamics 
depends on the sign of the chemical potential fi (which indicates how many particles 
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are trapped in the BEC) and is illustrated in figure |U Here, we considered the 
two situations with negative chemical potential (figure HJ),c). In both cases, there 
is a family of invariant tori winding around the trivial periodic orbit R = 0. When 
viewed in the proper coordinates, the dynamics on each torus is a flow with constant 
velocity vector. One frequency is equal to 1, and the other varies from one torus to 
the next, going to infinity monotonically as the distance of the torus to R = goes 
to infinity. The "proper coordinates" are the action-angle coordinates introduced 
in section 0] These tori correspond to quasiperiodic oscillations of R(x), which 
describes the amplitude dynamics of the BEC wave function. We note that in this 
unforced setting, the quasiperiodicity has no physical meaning, as the period of the 
optical lattice potential does not yet play a role. 

In the forced setting, V ^ 0. Figures [T] HH illustrate the case of a periodic lattice 
for a range of amplitudes. We remark that from a mathematical point of view, there 
is no extra difficulty in dealing with superlattices or even lattices with more than 
two Fourier modes as compared to lattices with just a single mode, provided the 
periods are commensurate (that is, provided they have a common multiple). Our 
main result, Theorem |3 uses KAM theory to demonstrate that for any size of the 
periodic forcing, there exist invariant tori winding around the trivial periodic orbit. 
Let -R m ax denote the maximal R coordinate on an invariant torus. That is, -R max 
is the amplitude of the corresponding oscillations. In terms of the mathematical 
setting Q, invariant tori exist for a dense set of "admissible" -R max in the interval 
(cz A l l 2 \ 1 22 1 1 1 / 2 , +oo), for some positive constant c > 0. In terms of the physical 
parameters, the lowest admissible amplitude is roughly equal to Sbihuoj / ^ — rag), 
where u is its frequency, k is the lattice wave number, b\ ~ 0.874019 is a constant, 
m is the mass of the atomic species in the condensate, and g < is the scaled value 
of the two-body scattering length. The parameter g varies from one condensate 
species to another and can be changed by exploiting Feshbach resonances [2T] . 
Authors' note: The following paragraph has been edited. 

The set of admissible -R max is characterized by a Diophantine condition on the 
frequencies of the torus. A slightly different KAM theorem than the one used in this 
work shows that, for each forcing, the relative measure of invariant tori converges 
to full measure as -R max — ► +oo. See, for example, |I7| EH] for a general statement 
or |H] for discussion of a system similar to ours. Thus, we find a large measure of 
quasiperiodic dynamics, and in the complement of the union of invariant tori, there 
exist Aubry-Mather sets. Furthermore, by the Poincare-Birkhoff Theorem [JJ E], 
in between any two invariant tori there are periodic orbits of saddle and center 
type for all intermediate resonant frequencies (technically, these can be classified as 
Aubry-Mather sets as well). As a result, there is a large measure of quasiperiodic 
invariant tori interlaced with Aubry-Mather sets and resonant layers. We further 
note that the homoclinic figure eight that exists in the unforced setting for u < 
will generically break under the forcing, generating a homoclinic tangle, as can be 
checked for this specific case by computing Melnikov integrals. We refer to Ell 
for such computations on similar systems. A possible direction for future study 
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would be to establish conditions for the nonexistence of invariant tori in the system 
under consideration, which would complement the existence conditions presented in 
this paper. Such a converse KAM theorem was obtained in [H] for area-preserving 
twist maps and in [39J in a Lagrangian/Hamiltonian setting. 

From a physical perspective, we recall that the GP equation is derived from 
a many-body quantum problem under the assumption that the Bose gas is dilute 
[TT| 13*3*] . In particular, this implies that the mean interparticle distance should be 
much larger than the scattering length. Therefore, the large R results obtained 
here would necessitate the condensate to have expanded sufficiently to ensure that 
its density is low (so that the reduction from the many-body problem to the three- 
dimensional GP description remains valid). Accordingly, for sufficiently large R, one 
would eventually have to include corrections to the GP equation arising beyond the 
mean-field description to describe the physics correctly. Incorporating beyond-mean- 
field dynamics in the study of BECs is a difficult problem (see, e.g., the discussion 
in Ref. |59j), and it is not agreed precisely when such corrections become relevant 
or what form they should take to simultaneously ensure tractability and correctly 
capture the physics. 
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